physica status solidi, 17 June 2010 



o 



Accuracy of Quantum Monte Carlo 
Methods for Point Defects in Solids 



William D. Parker^ John W. Wilkins\ Richard G. Hennig 



'1,2 



Department of Physics, The Ohio State University, 191 W. Woodruff Ave., Columbus, Ohio 43210, USA 
" Department of Materials Science and Engineering, Cornell University, Ithaca, New York 14853, USA 

Received XXXX, revised XXXX, accepted XXXX 
Published online XXXX 

PACS 



' Corresponding author: e-mail rhennig@cornell.eclu. Phone: +01-607-2546429 



Quantum Monte Carlo approaches such as the diffusion Monte Carlo (DMC) method are among the most accurate 
many-body methods for extended systems. Their scaling makes them well suited for defect calculations in solids. We 
review the various approximations needed for DMC calculations of solids and the results of previous DMC calculations 
for point defects in solids. Finally, we present estimates of how approximations affect the accuracy of calculations for 
self-interstitial formation energies in silicon and predict DMC values of 4.4(1), 5.1(1) and 4.7(1) eV for the X, T and 
H interstitial defects, respectively, in a 16(H-l)-atom supercell. 
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1 introduction Point defects, such as vacancies, in- 
terstitials and anti-site defects, are the only thermodynam- 
ically stable defects at finite temperatures |]JJ. The infinite 
slope of the entropy of mixing at infinitesimally small de- 
fect concentrations results in an infinite driving force for 
defect formation. As a result, at small defect concentra- 
tions, the entropy of mixing always overcomes the enthalpy 
of defect formations. In addition to being present in equi- 
librium, point defects often control the kinetics of mate- 
rials, such as diffusion and phase transformations and are 
important for materials processing. The presence of point 
defects in materials can fundamentally alter the electronic 
and mechanical properties of a material. This makes point 
defects technologically important for applications such as 
doping of semiconductors ||2l[3l, solid solution hardening 
of alloys |4,!5|, controlling the transition temperature for 
shape-memory alloys Q and the microstructural stabiliza- 
tion of two-phase superalloys. 

However, the properties of defects, such as their struc- 
tures and formation energies, are difficult to measure in 
some materials due to their small sizes, low concentra- 
tions, lack of suitable radioactive isotopes, etc. Quantum 
mechanical first-principles or ab initio theories make pre- 
dictions to fill in the gaps left by experiment ||7|- 

The most widely used method for the calculation of 
defect properties in solids is density functional theory 
(DFT). DPT replaces explicit many-body electron inter- 



actions with quasiparticles interacting via a mean-field 
potential, i.e. the exchange-correlation potential, which is 
a functional of the electron density ||8l|. A universally true 
exchange-correlation functional is unknown, and DFT cal- 
culations employ various approximate functionals, either 
based on a model system or an empirical fit. The most 
commonly used functionals are based on DMC simula- 
tions Q for the uniform electron gas at different densities, 
e.g., the local density approximation (LDA) IIIOIII II and 
gradient expansions, e.g., the generalized gradient approx- 
imation (GGA) III2III3lll4III51fT6l . These local and semi- 
local functionals suffer from a significant self-interaction 
error reflected in the variable accuracy of their predictions 
for defect formation energies, charge transition levels and 
band gaps II17II18II . Another class of functionals, called 
hybrid functionals, include a fraction of exact exchange to 
improve their accuracy III9II20I . 

The seemingly simple system of Si self-interstitials ex- 
emplifies the varied accuracy of different density func- 
tionals and many-body methods. The diffusion and ther- 
modynamics of silicon self-interstitial defects dominates 
the doping and subsequent annealing processes of crys- 
talline silicon for electronics applications 1 3, 2 1,22). The 
mechanism of self-diffusion in silicon is still under de- 
bate. Open questions Il23l include: (1) are the interstitial 
atoms the prime mediators of self-diffusion, (2) what is the 
specific mechanism by which the interstitials operate, and 
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(3) what is the value of the interstitial formation energy? 
Quantum mechanical methods are well suited to determine 
defect formation energies. LDA, GGA and hybrid func- 
tionals predict formation energies for these defects ranging 
from about 2 to 4.5 eV ll24ll . Quasiparticle methods such as 
the GW approximation reduce the self-interaction error in 
DFT and are expected to improve the accuracy of the inter- 
stitial formation energies. Recent Go Wo calculations ll25l 
predict formation energies of about 4.5 eV in close agree- 
ment with HSE hybrid functionals [1241 and previous DMC 
calculations li26ii24J . Quantum Monte Carlo methods pro- 
vide an alternative to DFT and a benchmark for defect for- 
mation energies II27II28I . 

In this paper, we review the approximations that are 
made in diffusion Monte Carlo (DMC) calculations for 
solids and estimate how these approximations affect the 
accuracy of point defect calculations, using the Si self- 
interstitial defects as an example. Section |2] describes the 
quantum Monte Carlo method and its approximations. Sec- 
tion|3]reviews previous quantum Monte Carlo calculations 
for defects in solids, and Section |4]discusses the results of 
our calculations for interstitials in siUcon and the accuracy 
of the various approximations. 

2 Quantum Monte Carlo method Quantum Monte 
Carlo (QMC) methods are among the most accurate elec- 
tronic structure methods available and, in principle, have 
the potential to outperform current computational meth- 
ods in both accuracy and cost for extended systems. QMC 
methods scale as 0{N'^) with system size and can han- 
dle large systems. At the present time, calculations for as 
many as 1,000 electrons on 1,000 processors make effec- 
tive use of available computational resources Il24ll . Current 
work is under way to develop algorithms that extend the 
system size accessible by QMC methods to petascale com- 
puters 1291 . 

Continuum electronic structure calculations primarily 
use two QMC methods Il27l : the simpler variational Monte 
Carlo (VMC) and the more sophisticated diffusion Monte 
Carlo (DMC). In VMC, a Monte Carlo method evaluates 
the many-dimensional integral to calculate quantum me- 
chanical expectation values. Accuracy of the results de- 
pends crucially on the quality of the trial wave function 
which is controlled by the functional form of the wave 
function and the optimization of the wave functions param- 
eters 130J. DMC removes most of the error in the trial wave 
function by stochastically projecting out the ground state 
using an integral form of the imaginary-time Schrodinger 
equation. 

One of the most accurate forms of trial wave functions 
for quantum Monte Carlo applications to problems in elec- 
tronic structure is a sum of Slater determinants of single- 
particle orbitals multiplied by a Jastrow factor and modi- 
fied by a backflow transformation: 

3'(r„) = e^^'--^'") ^c,;CSF,(a;„)- 



The Jastrow factor J typically consists of a low order 
polynomial and a plane-wave expansion in electron co- 
ordinates r„ and nuclear coordinates Rm that efficiently 
describe the dynamic correlations between electrons and 
nuclei. Static (near-degeneracy) correlations are described 
by a sum of Slater determinants. Symmetry-adapted lin- 
ear combinations of Slater determinants, so-called config- 
uration state functions (CSF), reduce the number of deter- 
minant parameters c,. For extended systems, the lack of 
size consistency for a finite sum of CSF's makes this form 
of trial wave functions impractical, and a single determi- 
nant is used instead. Finally, the backflow transformation 
Tn Xn allows the nodes of the trial wave function to 
be moved which can efficiently reduce the fixed-node er- 
ror lISTI . Since the backflow transformed coordinate of an 
electron x„ depends on the coordinates of all other elec- 
trons, the Sherman-Morrison formula used to efficiently 
update the Slater determinant does not apply, increasing 
the scaling of QMC to 0(7^). If a finite cutoff for the 
backflow transformation is used, the Sherman-Morrison- 
Woodbury formula 1321 applies and the scaling is reduced 
to 0{N^). 

Optimization of the many-body trial wave function is 
crucial because accurate trial wave functions reduce statis- 
tical and systematic errors in both VMC and DMC. Much 
effort has been spent on developing improved methods for 
optimizing many-body wave functions, and this continues 
to be the subject of ongoing research. Energy and vari- 
ance minimization methods can effectively optimize the 
wave function parameters in VMC calculations 130113311 . 
Recently developed energy optimization methods enable 
the efficient optimization of CSF coefficients and orbital 
parameters in addition to the Jastrow parameters for small 
molecular systems, eliminating the dependence of the re- 
sults on the input trial wave function l30l . 

VMC and DMC contain two categories of approxi- 
mation to make the many-electron solution tractable: con- 
trolled approximations, whose errors can be made arbitrar- 
ily small through adjustable parameters, and uncontrolled 
approximations, whose errors are unknown exactly. The 
controlled approximations include the finite DMC time 
step, the finite number of many-electron configurations that 
represent the DMC wave function, the basis set approxi- 
mation, e.g., spline or plane-wave representation, for the 
single-particle orbitals such of of the trial wave function 
and the finite-sized simulation cell. The uncontrolled ap- 
proximations include the fixed-node approximation which 
constraints the nodes of the wave function in DMC to be 
the same as the ones for the trial wave function, the re- 
placement of the core electrons around each atom with a 
pseudopotential to represent the core-valence electronic in- 
teraction and the locality approximation that uses the trial 
wave function to project the nonlocal angular momentum 
components of the pseudopotential. 

2.1 Controlled approximations 
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Time step Diffusion Monte Carlo is based on the 
transformation of the time-dependent Schrodinger equa- 
tion into an imaginary-time diffusion equation with a 
source-sink term. The propagation of the SA^-dimensional 
electron configurations (walkers) that sample the wave 
function requires a finite imaginary time step which intro- 
duces an error in the resulting energy 04113511 . 

Controlling the time step error is simply a matter of 
performing calculations for a range of time steps either to 
determine when the total energy or defect formation energy 
reaches the required accuracy or to perform an extrapola- 
tion to a zero time step using a low order polynomial fit of 
the energy as a function of time step. Smaller time steps, 
however, require a larger total number of steps to sample 
sufficiently the probability space. Thus, the optimal time 
step should be small enough to add no significant error to 
the average while large enough to keep the total number of 
Monte Carlo steps manageable. In addition, the more ac- 
curate the trial wave function is the smaller the error due to 
the time-step will be Il35ll . 

Configuration population In DMC, a finite number 
of electron configurations represent the many-body wave 
function. These configurations are the time-independent 
Schrodinger equation's analogues to particles in the dif- 
fusion equation and have also been called psips ll34l and 
walkers IZTl . To improve the efficiency of sampling the 
many-body wave function, the number of configurations is 
allowed to fluctuated from time step to time step in DMC 
using a branching step. However, the total number of con- 
figurations needs to be controlled to avoid the configura- 
tion population to diverge or vanish Il35l . This population 
control introduces a bias in the energy. In practice where 
tested 1361 . few hundreds of configurations are sufficient 
to reduce the population control bias in the DMC total en- 
ergy below the statistical uncertainty. 

The VMC and DMC calculations parallelize easily 
over walkers. After an initial decorrelation run, the prop- 
agation of a larger number of walkers is computationally 
equivalent to performing more time steps. The variance of 
the total energy scales like 

2 ^corr 
OC 

Nconl Nstep 

where A^conf denotes the number of walkers, iVstop the 
number of time steps and Tcorr the auto correlation time. 

Basis set A sum of basis functions with coefficients 
represents the single-particle orbitals in the Slater determi- 
nant. A DFT calculation usually determines these coeffi- 
cients. Plane waves provide a convenient basis for calcula- 
tions of extended systems since they form an orthogonal 
basis that systematically improves with increasing num- 
ber of plane waves that span the simulation cell. Increasing 
the number of plane waves until the total energy converges 
within an acceptable threshold in DFT creates a basis set 
that has presumably the same accuracy in QMC. 

Since the plane wave basis functions are extended 
throughout the simulation cell, the evaluation of an orbital 



at a given position requires a sum over all plane waves. 
Furthermore, the number of plane waves is proportional to 
the volume of the simulation cell. The computational cost 
of orbital evaluation can significantly be reduced by using 
a local basis, such as B-splines, which replaces the sum 
over plane waves with a sum over a small number of local 
basis functions. The resulting polynomial approximation 
reduces the computational cost of orbital evaluation at a 
single point from the number of plane waves (hundreds 
to thousands depending on the basis set) to the number 
of non-zero polynomials (64 for cubic splines) 1371 . The 
wavelength of the highest frequency plane wave sets the 
resolution of the splines. Thus, the most important quantity 
to control in the basis set approximation is the size of the 
basis set. 

Simulation cell Simulation cells with periodic bound- 
ary conditions are ideally suited to describe an infinite 
solid but result in undesirable finite-size errors that need 
correction. There are three types of finite-size errors. First, 
the single-particle finite-size error arises from the choice 
of a single fc-point in the single-particle Bloch orbitals of 
the trial wave function. Second, the many-body finite size 
error arises from the non-physical self-image interactions 
between electrons in neighboring cells. Third, the defect 
creates a strain field that results in an additional finite size 
error for small simulation cells. 

The single-particle finite size error is greatly reduced 
by averaging DMC calculations for single-particle orbitals 
at different fc-points that sample the first Brillouin zone of 
the simulation cell, so-called twist-averaging ||38]| Alterna- 
tively, the single-particle finite size error can also be esti- 
mated from the DFT energy difference between a calcu- 
lation with a dense A: -point mesh and one with the same 
single fc-point chosen for the orbitals of the QMC wave 
function. 

For the many-body finite size error, several methods 
aim to correct the fictitious periodic correlations between 
electrons in different simulation cells. The first approach, 
the model periodic Coulomb (MPC) interaction [|39l . re- 
vises the Ewald method |40| to account for the periodic- 
ity of the electrons by restoring the Coulomb interaction 
within the simulation cell and using the Ewald interaction 
to evaluate the Hartree energy. The second approach is 
based on the random phase approximation for long wave 
lengths. The resulting first-order, finite-size-correction 
term for both the kinetic and potential energies can be 
estimated from the electronic structure factor (41 )■ The 
third approach estimates the many-body finite-size er- 
ror from the energy difference between DFT calculations 
using a finite-sized and an infinite-sized model exchange- 
correlation functional ||42|| . This approach relies on the 
exchange-correlation functional being a reasonable de- 
scription of the system, whereas the other two approaches 
(MPC and structure factor) do not have this restriction. The 
MPC and structure factor corrections are fundamentally re- 
lated and often result in similar energy corrections 1431 . 
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The defect strain finite size error, can be estimated 
at the DFT level using extrapolations of large simulation 
cells. Also since QMC force calculations are expensive 
and still under development 1441 . QMC calculations for ex- 
tended systems typically start with DFT-relaxed structures. 
Energy changes due to small errors in the ionic position as 
well as thermal disorder are expected to be quite small be- 
cause of the quadratic nature of the minima and will largely 
cancel when taking energy differences for the defect ener- 
gies. 

2.2 Uncontrolled approximations 

Fixed-node approximation The Monte Carlo algo- 
rithm requires a probability distribution, which is non- 
negative everywhere, but fermions, such as electrons, are 
antisymmetric under exchange, and therefore any wave 
function of two or more fermions has regions of posi- 
tive and negative value. For quantum Monte Carlo to take 
the wave function as the probability distribution, Ander- 
son Il34l fixed the zeros or nodes of the wave function and 
took the absolute value of the wave function as the proba- 
bility distribution. If the trial wave function has the nodes 
of the ground state, then DMC projects out the ground 
state. However, if the nodes differ from the ground state, 
then DMC finds the closest ground state of the system 
within the inexact nodal surface imposed by the fixed-node 
condition. This inexact solution has an energy higher than 
that of the ground state. 

Three methods estimate the size of the fixed node ap- 
proximation: (1) In the Slater- Jastrow form of the wave 
function, the single-particle orbitals in the Slater deter- 
minant set the zeroes of the trial wave function. Since 
these orbitals come from DFT calculations, varying the 
exchange-correlation functional in DFT changes the trial 
wave function nodes and provides an estimate of the size 
of the fixed-node error. (2) Lopez Rios et al. OTI applied 
backflow to the nodes by modifying the interparticle dis- 
tances, enhancing electron-electron repulsion and electron- 
nucleus attraction. The expense of the method has thus far 
limited its application in the literature to studies of sec- 
ond and third-row atoms, the water dimer and the ID and 
2D electron gases. (3) Because the eigenfunction of the 
Hamiltonian has zero variance in DMC, a linear extrapo- 
lation from the variances of calculations with and without 
backflow to zero variance estimates the energy of the exact 
ground state of the Hamiltonian. 

Pseudopotential Valence electrons play the most 
significant roles in determining a composite system's prop- 
erties. The core electrons remain close to the nucleus and 
are largely inert. The separation of valence and core elec- 
tron energy scales allows the use of a pseudopotential to 
describe the core-valence interaction without explicitly 
simulating the core electrons. However, there is often no 
clear boundary between core and valence electrons, and the 
core-valence interaction is more complicated than a simple 
potential can describe. Nonetheless, the computational de- 
mands of explicitly simulating the core electrons and the 



practical success of calculations with pseudopotentials in 
reproducing experimental values promote their continued 
use in QMC. Nearly all solid-state and many molecular 
QMC calculations to date rely on pseudopotentials to re- 
duce the number of electrons and the time requirement of 
simulating the core-electron energy scales. 

Comparing DMC energies using pseudopotentials con- 
structed with different energy methods (DFT and Hartree- 
Fock[HF]) provides an estimate of the error incurred by 
the pseudopotential approximation. Additionally, the dif- 
ference between density functional pseudopotential and 
all-electron energies estimates the size of the error intro- 
duced by the pseudopotential and is used as a correction 
term. 

Pseudopotential locality DMC projects out the 
ground state of a trial wave function but does not produce 
a wave function, only a distribution of point-like config- 
urations. However, the pseudopotential contains separate 
potentials (or channels) for different angular-momenta of 
electrons. One channel, identified as local, does not require 
the wave function to evaluate, but the nonlocal channels 
require an angular integration to evaluate, and such an 
integration requires a wave function. Mitas et al. Il45l in- 
troduced use of the trial wave function to evaluate the 
nonlocal components requiring integration. This locality 
approximation has an error that varies in sign. While there 
are no good estimates of the magnitude of this error, Ca- 
sula |[46| developed a lattice-based technique that makes 
the total energy using a nonlocal potential an upper bound 
on the ground-state energy. Pozzo and Alfe [|47ll found that, 
in magnesium and magnesium hydride, the errors of the 
locality approximation and the lattice-regularized method 
are comparably small, but the lattice method requires a 
much smaller time step (0.05 Ha^^ vs. 1.00 Ha~^ in Mg 
and 0.01 Ha^^ vs. 0.05 Ha^^ in MgH2) to achieve the 
same energy. Thus, they chose the nonlocal approxima- 
tion. 

While all-electron calculations would, in principle, 
make the pseudopotential and locaUty errors controllable, 
in practice, the increase in number of electrons, required 
variational parameters and variance of the local energy 
makes such calculations currently impractical for anything 
but small systems and light elements 14811 . 

3 Review of previous DMC defect calculations 

To date, there have been DMC calculations for defects in 
three materials: the vacancy in diamond, the Schottky de- 
fect in MgO and the self-interstitials in Si. 

3.1 Diamond vacancy Diamond's high electron and 
hole mobility and its tolerance to high temperatures and ra- 
diation make it a technologically important semiconductor 
material. Diffusion in diamond is dominated by vacancy 
diffusion l53ll . and the vacancy is also associated with radi- 
ation damage l54l . Table[T]shows the range of vacancy for- 
mation and migration energies calculated by LDA l50l and 
DMC l49l. DMC used structures from LDA relaxation and 
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Figure 1 DMC, GW and DFT energies (in eV) for neutral defects in three materials. DMC and experimental values have 
an estimated uncertainty indicated by numbers in parenthesis. For the diamond vacancy, DFT-LDA and DMC include a 
0.36 eV Jahn-Teller relaxation energy. LDA relaxation produced the structures and transition path so the DMC value for 
migration energy is an upper bound on the true value. The Schottky energy in MgO is the energy to form a cation-anion 
vacancy pair DFT-LDA produces a range from 6-7 eV depending on the representation of the orbitals and treatment of the 
core electrons. DMC using a plane-wave basis and pseudopotentials results in a value on the upper end of the experimental 
range. For Si interstitial defects, DFT values of the formation energy range from 2 eV below up to the DMC values, 
depending on the exchange-correlation functional(LDA, GGA[PBE] or hybrid[HSE]), and the GW values lie within the 
two-standard-deviation confidence level of DMC. 
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single-particle orbitals employing a Gaussian basis. LDA 
pseudopotential removed core electrons. The DMC calcu- 
lations predict a lower formation energy than LDA. The 
DMC value for the migration energy is an upper bound 
on the actual number since the structures have not been 
relaxed in DMC. Furthermore, DMC estimates the experi- 
mentally observed dipole transition and provides an upper 
bound on the migration energy. ||49l The GRl optical tran- 
sition is not a transition between one-electron states but be- 
tween spin states and ^T2. DMC calculates a transition 
energy of L5(3) eV from ^E to ^T2, close to the experi- 
mentally observed value of 1.673 eV. LDA cannot distin- 
guish these states. For the cohesive energy, DMC predicts 
a value of 7.346(6) eV in excellent agreement with the ex- 
perimental result of 7.371(5) eV while LDA overbinds and 
yields 8.61 eV. 

3.2 MgO Schottky defect MgO is an important test 
material for understanding oxides. Its rock-salt crystal 
structure is simple, making it useful for computational 
study. Schottky defects are one of the main types of defects 
present after exposure to radiation, according to classical 
molecular dynamics simulations 1551 . Table [T] shows that 
DMC predicts a Schottky defect formation energy in MgO 
at the upper end of the range of experimental values ifSTI . 

3.3 Si interstitial defects Table [l] shows that DFT 
and DMC differ by up to 2 eV in their predictions of the 
formation energies of these defects i26..24 1 . We compare 
the DMC values with our results including tests on the 
QMC approximations in Section]?] 

4 Results We specifically test the time-step, pseu- 
dopotential and fixed-node approximations for the for- 
mation energies of three silicon self-interstitial defects, 
the split-(llO) interstitial (X), the tetrahedral interstitial 



(T) and the hexagonal interstitial (H). The QMC calcula- 
tions are performed using the CASINO 1561 code. Density 
functional calculations in this work used the QUANTUM 
ESPRESSO 157} and WIEN2k f5E\ codes. The defect 
structures are identical to those of Batista et al. [24|. 
The orbitals of the trial wave function come from DFT 
calculations using the LDA exchange-correlation func- 
tional. The plane-wave basis set with a cutoff energy of 
1 , 088 eV (60 Ha) converges the DFT total energies to 
1 meV. A 7x7x7 Monkhorst-Pack A: -point mesh centered 
at the L-point (0.5,0.5,0.5) converges the DFT total energy 
to 1 meV. A population of 1 , 280 walkers ensured that 
the error introduced by the population control is negligi- 
ble small. Due to the computational cost of backflow, we 
perform the simulations for a supercell of 16(h-1) atoms 
and estimate the finite-size corrections using the structure 
factor method j4li|. The final corrected DMC energies for 
the X, T and H defects are shown in the bottom line of 
Table]2] 

4.1 Time step Figure ]3] shows the total energies of 
bulk silicon and the X defect as a function of time step 
in DMC. A time step of 0.01 Ha~^ reduces the time step 
error to within the statistical uncertainty of the DMC total 
energy. 

4.2 Pseudopotential In our calculations, a Dirac- 
Fock (DF) pseudopotential represents the core electrons 
for each silicon atom I59II60I16T1 . To estimate the error 
introduced by the pseudopotential, we compare the defect 
formation energies in DFT using this pseudopotential with 
all-electron DFT calculations using the linearized aug- 
mented plane-wave method t58ll . This comparison gives 
corrections of 0.083, -0.168 and 0.054 eV for the H, T 
and X defects respectively. 
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Figure 2 DMC Si defect formation energies. Varying parameters and improved methods produce values for each defect 
that lie within two standard deviations of each other although the energetic ordering of the defects varies. All calculations 
use DFT-LDA to produce the orbitals in the Slater determinant. 
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Figure 3 DMC total energies with varying (imaginary) 
time steps for bulk silicon and the X defect. The error due 
to the finite-sized time step is smaller than the statistical 
uncertainty in the total energies for values from 10~^ to 
1 Ha^^. Note that these energies include no finite-size or 
pseudopotential corrections and thus differ in value from 
those in Figure |4] 



4.3 Fixed-node approximation The fixed-node ap- 
proximation is the main source of error of diffusion Monte- 
Carlo calculations. To estimate the size of the error intro- 
duced by the fixed-node approximation, we perform the 
calculations using the backflow transformation which al- 
lows the nodes of the trial wave function to be moved and 
reduces the fixed-node error |f3T|. We estimate the error 
due to the fixed-node approximation by performing calcu- 
lations with and without the backflow transformation and 
by extrapolating the resulting defect formation energies to 
zero variance. 



Defects-^ 



Bulk?- 




no backflow 



q-1814 

1815 
M816 

1817 
^-1818 

1819 
^-1820 



0.2 



0.4 



0.6 

'(eV) 



0.8 



1 



Figure 4 DMC total energies calculated with and without 
backflow transformation and linearly extrapolated to zero 
variance. The backflow transformation reduces the total en- 
ergies by 0.19(3), 0.62(5), 0.25(5) and 0.35(5) eV for bulk 
(solid line) and the X (dash-dot line), T (dashed line) and H 
(dotted line) defects respectively. The extrapolation to zero 
variance only slightly reduces the total energies by about 
0.1 eV. The reduction of the total energy from the back- 
flow transformation indicates the size of the errors due to 
the fixed-node and pseudopotential locahty approximation. 



Applying the backflow transformation to electron co- 
ordinates using polynomials of electron-electron, electron- 
nucleus and electron-electron-nucleus separations, we in- 
clude polynomial terms to eighth-order for each spin type 
in electron-electron separation, to sixth-order in electron- 
nucleus separation and to third-order for each spin type in 
electron-electron-nucleus separation. Figure |4] shows the 
linear extrapolation of the DMC energies for the Slater- 
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Jastrow and Slater- Jastrow-backflow trial wave function to 
zero variance. The total energy decrease for the bulk and 
interstitial cells due to the backflow transformation ranges 
from 0.20(5) to 0.62(5) eV. The backflow transformation 
results in a significantly improved nodal surface of the trial 
wave function which is reflected in the reduced the vari- 
ance of the local energy. 

Table |2] list the Si interstitials formation energies in 
DMC for the Slater-Jastrow, Slater-Jastrow-backflow wave 
function and the extrapolation. Applying the backflow 
transformation reduces the formation energies for the X, 
T and H interstitial by 0.42(5), 0.05(5) and 0.15(5) eV, 
respectively. The linear extrapolation provides a simple 
estimate of the remaining fixed node error The extrapola- 
tion lowers the interstitial formation energy by a negligible 
amount of 0.06(10), 0.00(10), 0.07(10) eV for the X, T 
and H interstitial, respectively. The resulting Si interstitial 
formation energies are 4.4(1), 5.1(1) and 4.7(1) eV for 
the X, T and H interstitial, respectively, in close agree- 
ment with recent GqWo ||25l and previous HSE and DMC 
calculations 112611241 . 

5 Conclusion QMC methods present an accurate tool 
for the calculation of point defect formation energies, pro- 
vided care is taken to control the accuracy of all the un- 
derlying approximations. Including corrections for the ap- 
proximations yields DMC values for the Si interstitial de- 
fects on par with GW and hybrid-functional DFT calcu- 
lations. While backflow transformation and zero-variance 
extrapolation to remove the fixed-node error modify the en- 
ergies slightly, further work remains to carefully control for 
finite-size effects known to plague defect supercell calcu- 
lations. 
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